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Abstract 

We discuss the irreversibility, nonlocality, and fluctuations, as well as the Lyapunov and hydro- 
dynamic instabilities characterizing atomistic, smooth-particle, and finite-difference solutions of 
the two-dimensional Rayleigh-Benard problem. To speed up the numerical analysis we control the 
time-dependence of the Rayleigh number, TZ(t) , so as to include many distinct flow morphologies 
in a single simulation. The relatively simple nature of these computational algorithms and the 
richness of the results they can yield make such studies and their interpretation particularly well 
suited to graduate-level research. 
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I. INTRODUCTION 



"Understanding Turbulence" is an enduring catch phrase and has been a potential funding 
source since the early days of computers. There is no shortage of reviews ranging from short 
sketches-^™ to scholarly studies^ - —. The vast research literature takes in earth, air, fire, and 
water as well as the weather, the sun, aircraft design, and small-box chaos. Spectra and 
power-law relations abound. Mostly the working fluid is incompressible and often its motion 
is described as a superposition of modes or vortices. Two- and three-dimensional turbulence 
behave differently, with the flow of energy away from or toward smaller length scales in these 
two cases^. "Enstrophy", the squared vorticity [ squared rotation rate ], is "conserved" in 
two-dimensional incompressible flow^. 

Despite all this information there appears to be more to learn. How many vortices should 
we expect to see? What is the Lyapunov spectrum like? How localized are the vectors 
corresponding to the exponents? The simple nature of the underlying model, a conducting 
viscous fluid, the complexity of the flows that result, and the multitude of computational 
schemes, all provide opportunities for imaginative approaches and analyses. We recommend 
their study and describe our own explorations of what seems to us the simplest problem 
involving turbulence, Rayleigh-Benard flow-i~— . 

The classical Rayleigh-Benard problem describes the convective behavior of a compress- 
ible, heat conducting, viscous fluid in the presence of gravity and a temperature gradient. 
Here we suppose that the fluid is confined by a stationary square L x L box with fixed 
boundary temperatures. Despite these simplest possible of boundary conditions, even in 
two space dimensions this problem provides interesting internal flows of mass, momentum, 
and energy. The heat driving these flows enters and exits along the boundaries. Most of 
it comes in at the bottom and flows out at the top. There are two competing mechanisms 
for the heat flow from bottom to top. The simpler of the two is conduction, described by 
Fourier's law, Q = — kVT . But mechanical ( or "convective" ) heat flow is possible too and 
comes to dominate conduction as the flow begins to move, and continues to grow as the flow 
eventually becomes turbulent. 

Thermal expansion near the bottom of the box provides the buoyancy necessary to carry 
the hot fluid upward. Cooling and compression near the top encourages downward flow. 
These vertical driving forces due to temperature and gravity are balanced by the dissipative 
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effects of heat conduction and viscosity which lead to macroscopic entropy production. The 
dimensionless ratios of these effects, the Rayleigh Number 71 and, to a lesser extent, the 
Prandtl Number V [ which we set equal to unity in our work here ] : 

K = g(d\nV/dT) P ATL 3 /(uD) ; V = (u/D) , 

control the overall flow. The Nusselt number M completes the list of dimensionless flow 
variables. It is an observable rather than an input. J\f is simply the ratio of the ( time- 
averaged, if necessary ) vertical heat flux to the prediction of Fourier's law : 

J\f = (LQ v /kAT) . 

With our thermostated sidewalls the definition of the Nusselt Number is somewhat arbitrary. 
Entropy production is a more appropriate measure of our flows' separation from equilibrium, 
though we will not discuss those interesting results here for lack of space. 

The convective flow patterns characterizing Rayleigh-Benard flow can be stationary, peri- 
odic in time, or chaotic. It is often possible to observe qualitatively different solutions - dif- 
ferent numbers of convective rolls for instance - for the same external boundary conditions^. 
And at very "high" Rayleigh numbers [ on the order of a half million or more ] , chaotic flows 
never repeat. Chaotic solutions describe at least two distinct regimes of turbulence^, called 
"soft" and "hard" , and distinguished by the form of their fluctuations, Gaussian or exponen- 
tial respectively™^. The time scales associated with eddy rotation vary from seconds in the 
laboratory to aeons inside the earth and sun. The richness of Rayleigh-Benard flow patterns, 
even or especially in two dimensions, together with their illustration of the fundamentals 
of fluid mechanics, instability theory, nonlinear dynamics, and irreversible thermodynamics, 
makes these problems an ideal introduction to the use of numerical methods in computa- 
tional fluid dynamics^^. 

We choose to study here the simplest possible constitutive model, an ideal gas , 

PV = NkT = E = Nme — (9 In V/dT) P = (1/T) ; (3/k) = ^ [ p\n(T '/ ' p)dxdy . 

Jo Jo 

For simplicity we set Boltzmann's constant k, the particle mass m, and the overall mass 
density p equal to unity. We choose the hot and cold temperatures equal to 1.5 and 0.5 so 
that AT ~ T. Finally, we choose the gravitational constant g so as to give a constant-density 
solution of the equation of motion in the quiescent purely conducting case : 

-(dP/dy) = pg = -p(kdT/dy) { g = (kAT/mL) ee (1/L) ->• p ~ 1 , a constant } . 
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Three Stick-Boundary Vortices 

Figure 1: An L x L square with an idealized roll using the "stick" boundary condition [ zero 
velocity at the wall ] : u x oc + sin(y) cos(a;/2) ; u y oc — sin(x) cos(y/2) ; — tt < x,y < +tt , 
is shown at the center. These data are also reflected, both to the right and to the left, to make 
an array with three rolls. This illustrates a handy initialization technique to use in the search for 
stable multiroll solutions. 

With these simplifications a square [ N = L x L ]-cell system with a Rayleigh number 71 
and a Prandtl Number V is achieved by choosing the two constitutive properties, kinematic 
viscosity v and thermal diffusivity D , to satisfy the two definitions : 

{y/D)=V\ K=(L/u)(L/D) . 

Unlike experimentalists we computational scientists are not limited to physical materials, 
dimensions, or boundary conditions. We have the undoubted luxury that our transport 
coefficients ( as well as the gravitational acceleration and even the box size ) can all be 
time dependent if we like. In the simulations reported here we typically use time-dependent 
transport coefficients, chosen so that the Rayleigh number increases or decreases [ to check 
for hysteresis ] linearly with time. In this way a whole range of Rayleigh numbers, with 
varying roll numbers, kinetic energies, and Lyapunov exponents ( if the increase is carried 
out sufficiently slowly ), can all be obtained with a single simulation. 

For sufficiently large values of the Rayleigh number ( 4960 or more for the static fixed- 
temperature boundary conditions used here^ ) one or more viscous conducting rolls form 
and evolve with time. With the convective heat flow directed upward, and on the average 
balanced by the gravitational forces acting downward, either stationary, periodic, or chaotic 
flows can be achieved. Figure 1 shows how a simple single vortex can be used to construct 
initial conditions with one or more vortices. 
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Since the 1980s nominally steady-state solutions for such flows have been computed with 
three distinct methods: microscopic molecular dynamics together with particle-based and 
grid-based macroscopic simulation methods^ - —. The Smooth- Particle Applied Mechan- 
ics Method ( SPAM ) offers a welcome bridge between the microscopic and macroscopic 
approaches^. In SPAM the dynamics of macroscopic particles is governed by motion equa- 
tions including the macroscopic irreversible constitutive laws. But the form of those motion 
equations mimics that of the microscopic motion equations. In both cases the accelera- 
tions are based on summed-up contributions from neighboring pairs of particles. SPAM 
calculations can also be thought of as a finite-difference algorithm on an irregular grid. 

The resulting macroscopic flow patterns exhibit interesting solution changes as the 
Rayleigh number increases above the critical value of 4960 . The positions of the rolls' 
centers can exhibit both periodic and chaotic motion. Lyapunov exponents characterize the 
growth of the instabilities leading to chaotic motion- 5 ^. For continuum simulations with 
thousands of degrees of freedom the simplest calculation of the instabilities involves only the 
largest exponent. There is evidence that the "spectrum" of Lyapunov exponents is roughly 
linear [ and with a negative sum, due to the dissipative nature of continuum flows ]*. 

It turns out that the relative stability of particular flows depends upon the initial condi- 
tions. No known variational principle ( like maximum entropy, or minimum entropy produc- 
tion ) predicts which of the several flows is stable 11 . The various "principles" based on energy 
or entropy can be evaluated for these simulations. Intercomparisons of the three simulation 
methods can shed light on the dissipation described by the Second Law of Thermodynamics 
and the differing time reversibilities of the microscopic and macroscopic techniques. 

In Section II we summarize the continuum physics of fluid flow problems: mass, momen- 
tum, and energy conservation are always required. Shear flows and heat flows can result. In 
Section III we outline three numerical solution techniques and display some typical results. 
In Section IV we present our conclusions and suggest research directions useful for students. 

II. CONTINUUM MECHANICS AND RAYLEIGH-BENARD FLOW 

A physical description of any continuum flow necessarily obeys the conservation laws for 
mass, momentum, and energy : 

p = — pV • u ; pit = — V ■ P + pg ; pe = —Vu : P — V ■ Q . 
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The simplest derivation of these three laws uses an Eulerian coordinate system fixed in 
space. The summed-up fluxes of each conserved quantity through the surfaces of each of the 
L 2 computational cells, plus the internal gravitational contributions give the time rates of 
change in the cells. In the Rayleigh-Benard problem the boundary source terms introduce 
and extract energy at the bottom, along the sides, and at the top, while the gravitational 
momentum density source pg acts throughout the volume . The fluid's constitutive proper- 
ties - the pressure tensor P and the heat flux vector Q - are computed from the local state 
variables { p,u,e } and their gradients. P and Q are the momentum and energy fluxes in a 
coordinate system "comoving" with the local velocity u(r, t) . 

It is easy to solve the continuum flow laws by converting them to sets of ordinary differen- 
tial equations. These latter equations incorporate the linear phenomenological constitutive 
relations pioneered by Newton and Fourier, expressing pressure in terms of the symmetrized 
velocity gradient, and the heat flux vector in terms of the temperature gradient : 

P = (P cq - AV • u)I - 77 [ Vu + W ] ; Q = -kVT . 

k is the thermal conductivity. We set the bulk viscosity equal to zero ( appropriate for an 
ideal gas ) by setting A + 77 = so that the pressure tensor has the following form : 

Pxx = p c q - v[ (du x /dx) - (du y /dy) ] ; 

P yy = P cq ~ v[ (du y /dy) - {du x /dx) ] ; 

p xy = -v[ (du y /dx) + (du x /dy) } . 

With these constitutive relations specified we have a well-posed continuum problem ready 
to solve. 

Figure 2 shows observed stationary roll patterns typical of a Rayleigh-Benard flow with 
a gravitational force acting in the negative y direction and a temperature gradient resulting 
in heat convection in the positive y direction. The temperature and velocity are fixed on 
the horizontal and vertical boundaries, just as in the idealized one- and three-roll flows of 
Figure 1 . Higher values of the Rayleigh number result in solutions that form with three or 
more rolls, periodic roll motions, and finally chaotic motions. See References 6 and 7 . 

Figure 3 shows four typical chaotic Rayleigh-Benard velocity plots. The Rayleigh Number 
here is 800 000 . A mesh of 160 x 160 cells and 161 x 161 nodes was used. 
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Figure 2: Stationary roll patterns observed at Rayleigh Numbers of 10,000 and 32,000 with 
a Prandtl number of unity. The velocities are taken from Eulerian grid-based solutions of the 
conservation laws with linear Newton- Fourier constitutive relations. 

III. NUMERICAL METHODS FOR RAYLEIGH-BENARD FLOWS 
A. Particle Methods: Nonequilibrium Molecular Dynamics and SPAM 

Nonequilibrium molecular dynamics is a straightforward but limited method for 
studying Rayleigh- Benard flows. Though the method is both simple and fundamental, atom- 
istic particle studies have several disadvantages : first, the equation of state can't be specified 
in advance ( only the interatomic force law is given in molecular dynamics ) ; second, the 
number of degrees of freedom required to simulate convective rolls is either thousands ( in 
two dimensions^ ) or millions ( in three dimensions^ ) ; third, the time step in molecular 
dynamics simulations is a fraction of the collision time rather than the considerably larger 
time [ dt = (dx/c) , where c is the sound velocity ] given by the continuum Courant condi- 
tion. Finally, even with these large particle numbers and small time steps, the fluctuations 
in the atomistic simulations are so large that snapshots of nominally steady flows show large 
deviations from time averages. 

Mareschal and Rapaport and their coworkers*^ have studied two- and three-dimensional 
molecular dynamics systems, relatively large at the time they were carried out ( with 5000 
and 3,507,170 particles respectively ). In both these cases time averages were required. The 
simulations did confirm that these time averages of the atomistic flows closely matched the 
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160 x 160 Cells, R = 800 000 




Figure 3: Four snapshots of a turbulent flow computed with 160 x 160 computational cells. 

corresponding stationary continuum simulations. We have carried out a few corroborating 
simulations. In these we used thermostated boundaries composed of particles tethered to 
fixed lattice sites. Rather than obeying conservative Newtonian mechanics, here the bottom 
row of "hot" and top row of "cold" boundary particles separately follow thermostated equa- 
tions of motion with their temperatures controlled by the Nose-Hoover friction coefficients 7 

( Chot j Ccold ) : 

{ f = (F({ r })/m) - Chotr }hot ; { f = (F({ r })/m) - Ccoiar }coid • 

Particle escapes can be prevented by using a strong repulsive boundary potential to reflect 
any particle venturing "outside" the box. Figure 4 compares a time-averaged exposure of a 
typical molecular dynamics run with 23,700 particles to the final snapshot from the same 
simulation. 

Smooth Particle Applied Mechanics ( SPAM 8 ) is a macroscopic particle alterna- 
tive to molecular dynamics. SPAM simulations are based on the continuum constitutive 
relations rather than atomistic interatomic forces. Hundreds of particles, rather than thou- 
sands, can generate rolls, the timestep is much larger, and individual snapshots do reproduce 
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Figure 4: Simulation of Rayleigh-Benard flow with molecular dynamics. A snapshot, using 

smooth-particle averages of the particle velocities, is at the left. Averages appear at the right. 

Normalized, Two Continuous Derivatives 



Lucy's Weight Function for < r < h 

Figure 5: Lucy's weight function, normalized for two dimensions, Jq 2irrw(r < h)dr = 1 . 

the stable roll structures quite well. SPAM defines local continuum averages by combin- 
ing contributions from a few dozen nearby particles. All of these continuum properties, 
{ p,u,e,P,Q, . . . } are local averages from sums using a weight function like Lucy's, which 
is shown in Figure 5 : 

w (r <h) = (5/7r/i 2 )(l + 3z)(l - zf ; z= (r/h) . 

Averages computed using this twice-differentiable weight function have two continuous spa- 
tial derivatives, enough for representing the righthand sides of the diffusive continuum equa- 
tions with continuous functions. 

Consider the simplest application of SPAM averaging, the definitions of the local densities 
and velocities in terms of smooth-particle weighted sums : 

p{r) = m ^ w(r — r,j) ; p(r)u(r) = m ^ w(r — rAvi . 

i i 

These definitions satisfy the continuum continuity equation exactly! The variation of the 
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density at a fixed location r can be evaluated by the chain rule : 

(dp/dt) r = m'^w'vi ■ (r% — r)/|r — . 

i 

Then notice that the gradient with respect to r of the product (pu) includes exactly the 
same terms, but with the opposite signs : 

V r ■ (pu) = m'^w'vi • (r — ri)/\r — ri\ . 

i 

Thus the SPAM version of the continuity equation , 

(dp/dt) = -V • (pu) < — > p = -pV ■ u , 

is an identity, independent of the form or range h of the weight function w(r < h) . This is 
not entirely a surprise as there is no ambiguity in the locations of the particles' masses and 
momenta. 

The pressure and energy are more complicated. The smooth-particle equation of motion^ 
is antisymmetric in the particle indices. Thus that motion equation , 

Vi = -m2 (P/P% + (P/P 2 ), ] ■ V iW (ri ~ r 3 ) ; 
j 

conserves linear momentum ( but not angular momentum ) exactly. Notice that whenever 
the pressure varies slowly in space the weight function plays the role of a repulsive potential 
with the strength of the interparticle "forces" proportional to the local pressure. 

The gradients in SPAM are evaluated by taking derivatives of the corresponding sums. 
The temperature gradient, for example, is : 

(VT)i = mJ2( T j ~ T iWijl ( r j ~ r i)/(\ r ij\Pij) } 5 Pij = VNpj or (p z + Pj)/2 . 
j 

Notice that two neighboring particles make no contribution to the temperature gradient if 
their temperatures match. With the gradients defined the pressure tensor and heat flux 
vectors can be evaluated for all the particles and used to advance the particle properties to 
the next time step : 

{ f,v,e } — ► { r,v,e } 

In all, the SPAM method averages involve about two dozen distinct particle properties. This 
computational effort is compensated by SPAM's longer length and time scales. 
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In addition to providing an alternative approach to solving the continuum equations the 
SPAM averaging technique can be used to average molecular dynamics properties such as P 
and Q. This approach is particularly valuable in Shockwaves ( see Chapter 6 of Reference 7 ), 
where constitutive properties change on an atomistic distance scale. It would be interesting 
to compare the two sides of the continuum energy and motion equations and to use this 
comparison to optimize the choice of the weight function's range h . 

The molecular dynamics results differ qualitatively from continuum results in their time 
symmetry, so that averaging offers a way of reducing this conflict. Time-dependent solutions 
offer a specially flexible technique for bringing the two approaches into better agreement. 
The weight functions also offer a way of carrying out the coarse graining which could be 
used to reduce the conflict between the microscopic and macroscopic forms of mechanics. 

B. Eulerian Finite-Difference Method 9 

Straightforward centered-difference approximations to the continuum equations provide a 
useful approach to the Rayleigh-Benard problem. Mareschal and his coworkers pointed out 
that an efficient numerical method can be based on square cells or zones, with the velocities 
and energies defined at the nodes and the densities defined in the cells. A small 10 x 10 cell 
program written in this way would solve 3x11x11 + 10 x 10 = 463 ordinary differential 
equations. The solution procedure follows a seven-step plan: [1] use linear interpolation and 
extrapolation near the boundaries to find the complete set of 484 nodal variables and 400 cell 
variables; [2] use centered differences to find Vm and VT ; [3] use these gradients to obtain 
P and Q ; [4] evaluate Vm : P and V • Q ; [5] evaluate (dp/dt) from the neighboring nodal 
values ; [6] evaluate (du/dt) from the pressure gradients and (de/dt) by summing the con- 
vective contributions and the work and heat; [7] use fourth-order Runge-Kutta integration 
to advance the 463 dependent variables to the next time step. 

The Rayleigh-Benard solutions - simple rolls 1 -, periodic solutions, or chaos - can be 
observed in either two dimensions, where there are plenty of puzzles to solve, or three. 
Because simulation and visualization are simpler in two dimensions, while the challenges to 
understanding remain severe, we choose two dimensions. The critical Rayleigh Number of 
about 5000 corresponds to an eddy width which can easily be resolved with 8x8 cells and 
a one-roll Reynolds number of order unity. 
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Roll -> Rolls -> Periodic --> Chaotic 




< R = time < 700 000 

Figure 6: Kinetic energy per cell (vertical at top and horizontal at bottom) from a simulation 
with TZ = t. Notice, at the extreme left of the two plots, the reduction in the horizontal kinetic 
energy at the transition from one roll to two ( TZ ~ 25K ) . 

The Rayleigh number varies as L 4 . Doubling the sidelength L — > 2L with g and the 
transport coefficients fixed changes TZ : 5000 — > 80000 and increases the number of rolls to 
four. Desktop or laptop machines are quite capable of simulations with 1Z = 1 000 000 , for 
which this simple-minded reasoning could lead us to expect about (1000000/5000) 1 / 4 ~ 14 
rolls. In fact this doesn't happen. See Figure 7 below. In two dimensions the energy flow 
is toward, rather than away from, large rolls. For 7Z = 800i^ [ K indicates thousands ] and 
a 160 x 160 mesh one finds occcasional deep minima in the time- dependent kinetic energy. 
These minima correspond to only two large rolls, as in the simple solutions without chaos, 
with TZ ~ lOfT . In three dimensions the chaotic flow is qualitatively different, and more 
complicated. Instead of whirling vortices one finds plumes ascending and descending, with 
mushroom shaped heads for large Prandtl numbers ( glycerin ) where viscosity dominates 
conductivity^. 

Figure 6 shows the variation of the kinetic energy per cell with the Rayleigh Number, 
where TZ = t < 700K. The transitions go from one roll to two, and from two rolls to a 
time-periodic arrangement with perhaps four, which in turn gives rise to chaos. Our simple 
centered-finite-difference fixed-timestep code "blew up" at Rayleigh numbers of 

{ 715K, 810JT, 8A0K, 860K, 905fsT, 9A0K } for L = { 16, 24, 32, 48, 64, 96 } . 

Figure 6 suggests that chaos sets in around TZ = 385K and gradually increases in strength 
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Six 24 x 24 Cell Ideal-Gas Snapshots 




Figure 7: Snapshots at time 0, 100, 200 (bottom row, left to right) and 300, 400, 500 (top row, 
left to right) starting with an instantaneously reversed chaotic flow with 1Z = 800 000 . 

until the algorithm becomes unstable for the chosen mesh. 

The motion responds relatively quickly to perturbations. To demonstrate this we show in 
Figure 7 six snapshots from a 24 x 24 simulation with all of the flow velocities instantaneously 
reversed from forward to backward at time 0, where the forward chaotic flow has three 
distinct rolls. By a time of 500 (a few dozen sound traversal times) the flow is back to 
normal, again with three rolls. But in the transformational process of discarding the time- 
reversed morphology the stabilizing flow acquires as many as seven rolls, in accord with the 
estimate that roughly 8x8 cells are required to resolve a roll. The quick reduction in roll 
number as steady chaos is regained is again a symptom of the flow of energy from smaller 
to larger rolls in two dimensions. 

Although there is no difficulty in simulating the motion on larger meshes even the modest 
576-cell simulation of Figure 7 provides an excellent characterization of chaos. Character- 
izing chaos quantitatively entails evaluating Lyapunov exponents, the tendency of nearby 
trajectories to separate farther or to approach one another. The separation of a satellite 
trajectory from its reference has four component types : 

$ = ( {M> {Su x }, {Su y }, {5e} ) . 
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Figure 8: Separation between two strongly chaotic simulations. The fourth-order Runge-Kutta 
timestep is 0.05. Although four different sets of data (for p, u, e subspaces as well as the complete 
space) are shown, here, the local Lyapunov exponents are scarcely distinguishable from one another. 

Keeping the separation constant by rescaling at every timestep gives the local exponent , 

\i(t) = ln(| <5 before | /| <5 aftcr \)/dt . 

In atomistic simulations it is well known that the largest Lyapunov exponent can be 
measured in configuration, momentum, or phase space, with identical longtime averages^. 
In continuum simulations one would expect that the density, velocity, or energy subspaces 
could be used in this same way. By choosing particular mass, length, and time units any 
one of the three subspaces could be made to dominate the local Lyapunov exponent. 

Figure 8 shows the exponential separation of a satellite trajectory from its reference 
trajectory, without rescaling, as measured in the density, velocity, energy, and complete state 
spaces of the flow. The similarity of the four curves is so complete that we don't attempt 
to label them separatedly in the figure. Evidently, despite the changing morphology of the 
chaotic vortices, the underlying chaos ( at a Rayleigh Number of 800K ) is quite steady. 
The research literature indicates that the Lyapunov spectrum in such two-dimensional flows 
is roughly linear, and corresponds to a strange attractor with only a few degrees of freedom, 
no doubt corresponding to the number of observed vortices. 

The low dimensionality of two-dimensional turbulent chaos results from the tremendous 
dissipation inherent in the Navier- Stokes-Newton- Fourier model. If we evaluate the three 
"phase-space derivatives" which contribute to the continuum analog of Liouville's particle 
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Theorem : 

(dp I dp) ; (dii/du) ; (de/de) , 

[ here the dots are time derivatives at fixed cells or nodes ] the velocity and energy derivatives 
give —4(r)/p)/(dx) 2 and — 4k/ (dx) 2 respectively while the density derivative vanishes. In the 
end only a few degrees of freedom exhibit chaos. 

IV. CONCLUSIONS AND SUGGESTIONS FOR RESEARCH 

Continuum mechanics can be studied with finite-difference ordinary differential equations, 
or with particle differential equations, resembling those used in molecular dynamics. The 
finite-difference approach is certainly the most efficient of these possibilities. The Rayleigh- 
Benard problems exhibit a variety of flows, with interesting results at the level of a few 
hundred compuational cells. The relative stability of the flows and the characterization 
of the chaos are both interesting research areas. With the limited dimensionality of the 
chaotic flows' attractors, estimating only a few Lyapunov exponents suffices to characterize 
Rayleigh-Benard chaos. The loci of Lyapunov vectors' instability is yet another source of 
fascinating questions and answers. 

Though we have no space to discuss them here the useful smooth-particle technique 
for bridging together the particle and continuum methods suggests a variety of problems 
designed to reduce the discrepancies between the three types of numerical algorithm. 
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